library(readr)
library(dplyr)
library(leaflet)
library(sf)
library(tidyr)
library(purrr)
library(tibble)
library(RColorBrewer)
library(gridExtra)
library(geosphere)
library(lubridate)
library(patchwork)

1 Carga de datos

La idea de este punto es obtener para cada zona de contaminación, las medias de intensidad, ocupacion, velocidad… de cada uno de los tres tipos de vias que se encuentran en un radio dado

trafico_final <- read_csv('/Users/pablogandia/Desktop/trafico_cluster.csv')
contaminantes <- read_csv("/Users/pablogandia/Desktop/Analisis_aire/data/raw/geoespacial/estaciones_valencia.csv") %>%
  select(Estación, Latitud, Longitud, Altitud)

zonas <- c("València - Av. França",
           "València - Molí del Sol", "València - Pista de Silla",
           "València - Politècnic", "València Port llit antic Túria",
           "València Port Moll Trans. Ponent")

contaminantes <- contaminantes[contaminantes$Estación %in% zonas, ]

datos_vias <- read_csv('/Users/pablogandia/Desktop/trafico_horario_clean.csv') %>%
  select(IdPM, Hora, Dia,Mes,IntensidadMediana,VelocidadMedia,OcupacionMedia,OcupacionMaxima,DevIntensidad)

Aquí establecemos para visualizar un radio de 1km, que puede ser interesante para analizar relaciones cercanas pero tampoco obviar las vias que tiene a una distancia media, que pueden ser muchas:

pal <- colorFactor(palette = brewer.pal(3, "Set1"),
                   domain = trafico_final$tipo_via)

m <- leaflet(trafico_final) %>%
  addTiles() %>%
  setView(lng = mean(trafico_final$Lon), lat = mean(trafico_final$Lat), zoom = 11) %>%
  addCircleMarkers(
    lng = ~Lon, lat = ~Lat,
    radius = 6,
    color = ~pal(tipo_via),
    fillOpacity = 0.8,
    stroke = FALSE,
    popup = ~paste("Zona:", idpm, "<br>Tipo de vía:", tipo_via)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~tipo_via,
    title = "Tipos de vía",
    opacity = 1
  ) %>%
  addCircles(
    data = contaminantes,
    lng = ~Longitud,
    lat = ~Latitud,
    radius = 1000,
    color = "#ff9999",
    fillOpacity = 0.2,
    stroke = TRUE
  ) %>%
  addCircleMarkers(
    data = contaminantes,
    lng = ~Longitud,
    lat = ~Latitud,
    radius = 8,
    color = "black",
    fillOpacity = 1,
    popup = ~Estación
  )

m

Vemos que principalmente las 4 estaciones que mas nos interesan son El Politécncio, La Avenida de Francia, la Pista de Silla y Moli del sol.

Las estaciones del puerto, aunque en Valencia, se ecuentran bastante distanciadas de las zonas de tráfica, y además están sesgadas por el tráfico marítimo, lo que puede complicar los analisis.

gdf_vias <- st_as_sf(trafico_final, coords = c("Lon", "Lat"), crs = 4326)
gdf_cont <- st_as_sf(contaminantes, coords = c("Longitud", "Latitud"), crs = 4326)

gdf_vias_m <- st_transform(gdf_vias, 3857)
gdf_cont_m <- st_transform(gdf_cont, 3857)

radio <- 1000
gdf_cont_buf <- st_buffer(gdf_cont_m, dist = radio)

estacion_a_vias <- map2(
  gdf_cont_buf$geometry, gdf_cont$Estación,
  ~ {
    vias_dentro <- gdf_vias_m[st_within(gdf_vias_m, .x, sparse = FALSE), ]
    setNames(list(vias_dentro$idpm), .y)
  }
) %>% reduce(c)

estacion_a_vias[["València - Molí del Sol"]]
##  [1] 1702 1715 1716 1717 1718 1719 3601 3612 3613 3622 3624 3625 3638 3639 3640
## [16] 3822
vias_info <- trafico_final %>%
  select(idpm, tipo_via)

df_estacion_vias <- tibble(
  estacion = rep(names(estacion_a_vias), lengths(estacion_a_vias)),
  idpm     = unlist(estacion_a_vias)
)

df_estacion_vias <- df_estacion_vias %>%
  left_join(vias_info, by = "idpm") %>%
  mutate(tipo_via = factor(tipo_via, levels = c("alto", "bajo", "congestionado")))

contadores <- df_estacion_vias %>%
  count(estacion, tipo_via) %>%
  pivot_wider(
    names_from = tipo_via,
    values_from = n,
    values_fill = 0,
    names_prefix = "tipo_"
  )

contadores

Aquí vemos la cantidad de vias de diferente tipo que se le asignan a cada estacion poniendo 1km de radio.

trafico <- trafico_final %>%
  left_join(datos_vias, by = c("idpm" = "IdPM"))

trafico <- trafico%>%
  select(-nombre,-Lat,-Lon)

trafico

2 Cruzar con contaminacion segun Estacion y tipo de vida

contaminacion <- read_csv('/Users/pablogandia/Downloads/Calidad_Aire_limpio.csv')

contaminacion <- contaminacion %>%
  mutate(
    Año = year(FechaHora),
    Mes = month(FechaHora),
    Dia = day(FechaHora),
    Hora = hour(FechaHora)
  )

zonas_interes <- c("València - Av. França",
                   "València - Molí del Sol",
                   "València - Pista de Silla",
                   "València - Politècnic",
                   "València Port llit antic Túria",
                   "València Port Moll Trans. Ponent")

contaminacion<- contaminacion %>%
  filter(Origen %in% zonas_interes, Año == 2023)
generar_base_estacion <- function(nombre_estacion, trafico, contaminacion, estacion_a_vias) {

  ids <- estacion_a_vias[[nombre_estacion]]
  
  trafico_filtrado <- trafico %>%
    filter(idpm %in% ids)
  
  global <- trafico_filtrado %>%
    group_by(Mes, Dia, Hora) %>%
    summarise(
      IntensidadGlobal = mean(IntensidadMediana, na.rm = TRUE),
      VelocidadGlobal = mean(VelocidadMedia, na.rm = TRUE),
      OcupacionGlobal = mean(OcupacionMaxima, na.rm = TRUE),
      .groups = "drop"
    )
  
  por_tipo <- trafico_filtrado %>%
    group_by(Mes, Dia, Hora, tipo_via) %>%
    summarise(
      Intensidad = mean(IntensidadMediana, na.rm = TRUE),
      Velocidad = mean(VelocidadMedia, na.rm = TRUE),
      Ocupacion = mean(OcupacionMaxima, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_wider(
      names_from = tipo_via,
      values_from = c(Intensidad, Velocidad, Ocupacion),
      names_glue = "{.value}_{tipo_via}",
      values_fill = 0
    )
  
  trafico_final <- left_join(global, por_tipo, by = c("Mes", "Dia", "Hora"))
  
  cont_filtrada <- contaminacion %>%
    filter(Origen == nombre_estacion) %>%
    mutate(
      Mes = month(FechaHora),
      Dia = day(FechaHora),
      Hora = hour(FechaHora)
    )

  base_final <- inner_join(cont_filtrada, trafico_final, by = c("Mes", "Dia", "Hora"))

  return(base_final)
}



base_franca <- generar_base_estacion("València - Av. França", trafico, contaminacion, estacion_a_vias)
base_franca

Así quedarían los datos cruzados con contaminación (misma zona y misma fecha), y ahora además le vamos a añadir la dirección del viento de los datos del objetivo 2.

colSums(is.na(base_franca))
##                FechaHora                   Origen                       NO 
##                        0                        0                        0 
##                      NO2                      NOx                       O3 
##                        0                        0                        0 
##                     PM10                    PM2.5                      SO2 
##                        0                        0                        0 
##                      Año                      Mes                      Dia 
##                        0                        0                        0 
##                     Hora         IntensidadGlobal          VelocidadGlobal 
##                        0                        0                        0 
##          OcupacionGlobal          Intensidad_alto          Intensidad_bajo 
##                        0                        0                        0 
## Intensidad_congestionado           Velocidad_alto           Velocidad_bajo 
##                        0                        0                        0 
##  Velocidad_congestionado           Ocupacion_alto           Ocupacion_bajo 
##                        0                        0                        0 
##  Ocupacion_congestionado 
##                        0

Algunos valores de cong y alta tienen nulos ya que se dan en zonas en las que no hay zonas de alta intensidad ni zonas de congestion frecuente. Por lo que estos nulos los imputaremos por 0.

3 Añadir direccion del Viento

viento <- read_csv('/Users/pablogandia/Downloads/meteorología_horaria_limpio.csv')

viento <- viento %>%
  select(fecha,hora,"Dirección del viento")

viento <- viento %>%
  mutate(
    fecha = as.Date(fecha),
    Mes = month(fecha),
    Dia = day(fecha),
    Hora = as.integer(hora)
  )


simplificar_direccion <- function(dir) {
  case_when(
    dir %in% c("N", "NNE", "NNW") ~ "N",
    dir %in% c("NE", "ENE") ~ "NE",
    dir %in% c("E", "ESE") ~ "E",
    dir %in% c("SE", "SSE") ~ "SE",
    dir %in% c("S", "SSW") ~ "S",
    dir %in% c("SW", "WSW") ~ "SW",
    dir %in% c("W", "WNW") ~ "W",
    dir == "NW" ~ "NW",
    TRUE ~ NA_character_
  )
}

viento <- viento %>%
  mutate(direccion = simplificar_direccion(`Dirección del viento`))

viento <- viento %>%
  select(Mes,Dia,Hora,direccion)

viento
base_franca <- base_franca %>%
  left_join(viento, by = c("Mes", "Dia", "Hora")) 
base_franca <- base_franca %>%
  mutate(
    viento_Norte     = as.integer(direccion %in% c("N", "NNE", "NNW")),
    viento_Sur       = as.integer(direccion %in% c("S", "SSW", "SSE")),
    viento_Este      = as.integer(direccion %in% c("E", "ESE", "ENE")),
    viento_Oeste     = as.integer(direccion %in% c("W", "WSW", "WNW")),
    viento_Noreste   = as.integer(direccion == "NE"),
    viento_Sureste   = as.integer(direccion == "SE"),
    viento_Suroeste  = as.integer(direccion == "SW"),
    viento_Noroeste  = as.integer(direccion == "NW"),
    viento_Sin_dato  = as.integer(is.na(direccion))
  )


base_franca

El resultado final es una base de datos que estudia contaminación en una zona y fecha concreta, asi como la direccion del viento y el tráfico de sus alrededores bastante bien carcterizado y segmentado.